Nonstationary spectrum estimation and time-frequency concentration

This paper extends Thomson's multitaper spectrum estimation method [17] to nonstationary signals. The method uses a newly{derived set of basis functions which generalize the concentration properties of the prolate spheroidal waveforms [15] to the time{frequency case. We solve for the basis which diagonalizes the nonstationary spectrum generating operator over a nite region of the time{frequency plane. These eigenfunctions are maximally concentrated to and orthogonal over the speci ed time-frequency region, and are thus doubly orthogonal. Individual spectrograms computed with these eigenfunctions form direct time{frequency spectrum estimates. We next present a multitaper time{ frequency spectrum estimation procedure using these time{ frequency eigenestimates. Bias and variance expressions are derived, allowing for a statistical characterization of the accuracy of the estimate. The time{frequency concentration property of the basis functions yields an estimator with excellent bias properties, while the variance of the estimate is reduced through the use of multiple orthogonal windows. 1. TIME{FREQUENCY SPECTRAL ANALYSIS There have generally been two approaches to time{frequency spectral analysis. The evolutionary spectrum approaches (e.g., [14, 7, 8]) model the spectrum as a slowly varying envelope of a complex sinusoid. This assumption allows the averaging of short{time spectral estimates to stabilize the variance. The second approach is commonly referred to as Cohen's bilinear class [3], which provides a general formulation for joint time{frequency distributions. Computationally, the evolutionary spectrum methods fall within Cohen's class. A subclass of time{frequency distributions are the positive time{frequency distributions (TFDs) [4]. Postive TFDs are everywhere nonnegative, and yield the correct univariate marginal distributions in time and frequency (the instantaneous energy and the energy spectral density): P (t; !) 0; (1) Z P (t; !)d! = js(t)j; (2) Z P (t; !)dt = jS(!)j; (3) where S(!) denotes the Fourier transform of the nite energy signal s(t), and all integrals are from 1 to 1. The rst method for generating positive TFDs used constrained optimization, minimizing the cross{entropy to a prior distribution subject to a set of linear constraints [9]. Positive TFDs have been linked to the evolutionary spectrum and estimated via deconvolution [13]. Least-squares estimation has also been used to compute positive TFDs [11]. Approximate solutions for positive TFDs have been obtained through a nonlinear combination of spectrograms [10]. Another approach to computing time{frequency spectra has been to extend Thomson's multitaper spectral estimation method [17] to the nonstationary case through a sliding{window framework [16]. [1] developed a multitaper time{frequency spectrum, including a signi cance test for nonstationary tones, using Hermite windows, which have previously been shown to maximize a time{frequency concentration measure [5]. [2] extended the Hermite multiwindow method to include a means of reducing artifacts using a time{frequency mask. While these methods all provide some representation of the time{varying frequency content of a signal, they do not relate the computed distribution to an underlying time{ frequency spectrum (e.g., [1] minimizes the bias between the multitaper TFD and the Wigner distribution; however, the Wigner distribution is not nonnegative for arbitrary signals, and as such is not a valid time{frequency spectrum). As a result, there is no quantitative measurement of the accuracy of the representation. For time{frequency analysis to be useful in a wide variety of real{world applications, some method of measuring the bias and variance of the estimated time{frequency spectrum is required. To meet this requirement, we present a statistical spectral estimation method for nonstationary signals. The method is based on a time{ varying lter formulation for positive TFDs, as discussed in [12]. We solve for the eigenvectors which diagonalize the nonstationary spectral generating function. These eigenvectors are maximally concentrated (and doubly{orthogonal) 1Throughout the analysis that follows, we use integral formulations of the various operations. The extension to the discrete, nite case is straightforward and not presented here. Integrals with no limits are over the entire domain of support of the integrand. The corresponding summations in the discrete case are then over the length of the corresponding vectors. in time{frequency. We then derive a multitaper estimation procedure to solve for the time{frequency spectrum. We also present bias and variance measures for the estimated time{frequency spectrum. 2. INTEGRAL EQUATION FOR A TIME{FREQUENCY SPECTRUM As is the case in stationary spectral estimation, a rigorous approach to time{frequency spectral estimation should be based upon the integral formulation underlying the generation of nonstationary signals. The formulation used here is a straightforward extension of the spectral representation theorem for stationary processes [14], and is equivalent to a linear time{varying (LTV) lter model. De ne the signal s(t) as the output of a white-noise-driven LTV lter: s(t) = Z h(t; )e( )d : (4) e(t) is bandlimited Gaussian white noise with bandwidth much greater than that of the lter h(t; ): e(t) = Z edZ(!): (5) dZ(!) is an orthogonal process with unit variance. The signal can then be written as: s(t) = Z H(t; !)edZ(!); (6) whereH(t; !) is de ned as the Fourier transform of h(t; t ) [12]. The time{frequency spectrum is de ned by: P (t; !) = jH(t; !)j: (7) This formulation for a time{frequency spectrum is of the same general form as Priestley's evolutionary spectrum [14]. However, we do not require that H(t; !) be slowly{varying. This form for P (t; !) also satis es the stochastic equivalent of the time and frequency marginals (equations 2{3); the relationship between the above time{varying spectrum and positive TFDs is discussed in [12]. Given a signal s(t), we want to estimate P (t; !); however, direct inversion of equation 6 is impossible. We can gain some idea of the time{varying frequency content of s(t) by computing the short-time Fourier transform (STFT): Ss(t; !) = Z s( )g(t )e j! d ; (8) where g(t) is a rectangular window of length T . The relationship between the STFT and H(t; !) is obtained by replacing s(t) by its time{frequency spectral formulation: Ss(t; !) = Z Z H( ; )g(t )e j(! ) dZ( )d : (9) To solve for the time-varying spectrum H( ; ), we need to invert the STFT operator g(t )e j! . This inversion is an inherently ill-posed problem. Instead, we approximate the inverse solution by regularizing it to some region R(t; !) in the time{frequency plane, much as Thomson regularized the spectral inversion to a bandwidth W in his multitaper approach [17]. For simplicity throughout, we will de ne R(t; !) to be a square region of time{frequency of dimension T W ; however, the results readily generalize to arbitrary regions R(t; !). In the case of spectral estimation, the operator is square and Toeplitz; its regularized inverse is found through an eigenvector decomposition. Such is not the case in the time{ frequency problem; the STFT operator is neither full rank nor square. To diagonalize it, we apply a Singular Value Decomposition, nding the left and right eigenvectors u( ) and V (t; !) and the associated eigen (singular) values :